import numpy as np
import matplotlib.pyplot as plt

# Set up the plot configurations.
plt.rcParams['lines.linewidth'] = 2.0
plt.rcParams['lines.color'] = 'black'
plt.rcParams['legend.frameon'] = True
plt.rcParams['figure.figsize'] = (8, 6)
plt.rcParams['font.family'] = 'serif'
plt.rcParams['legend.fontsize'] = 15
plt.rcParams['font.size'] = 15
plt.rcParams['axes.spines.right'] = False
plt.rcParams['axes.spines.top'] = False
plt.rcParams['axes.spines.left'] = True
plt.rcParams['axes.spines.bottom'] = True
plt.rcParams['axes.axisbelow'] = True
plt.rcParams['grid.color'] = 'grey'
plt.rcParams['grid.alpha'] = 0.6
plt.rcParams['grid.linestyle'] = '--'

sizex = 7.5
sizey = 5

# Create a figure and axis.
fig, ax = plt.subplots(2, 1,figsize=(sizex,2*sizey))

F=1
E=1
A=1
l=20
n=1

x = np.linspace(0, l, 100)


# Plotting the functions.
ax[0].plot(x, (F/(E*A)+n*l/E)*x-(n/(2*E))*x**2,color='red',linewidth=2.0,label="analytische Lösung")

fem_sol=np.array([
 [0,0],
 [l/2,F*l/(2*E*A)+3*n*l**2/(8*E)],
 [l,F*l/(E*A)+n*l**2 / (2*E)]   
])

#ax[0].scatter(fem_sol[:,0],fem_sol[:,1],color="blue",marker="o",s=100,label="FEM")
ax[0].plot(fem_sol[:,0],fem_sol[:,1],color="blue",marker="o",markersize=5,label="FEM")
# Set labels for the axis
ax[0].set_xlabel(r"$x$")
ax[0].set_ylabel(r"$u(x)$")
ax[0].set_title("Verschiebungsverlauf")

# Set custom ticks for the x-axis with LaTeX symbols
tick_positions = [0, l/2, l]
tick_labels = [r"$0$", r"$\frac{\ell}{2}$", r"$\ell$"]
ax[0].set_xticks(tick_positions)
ax[0].set_xticklabels(tick_labels)
ax[0].set_yticks([0,F*l/(2*E*A)+3*n*l**2/(8*E),F*l/(E*A)+n*l**2 / (2*E)])
ax[0].set_yticklabels([r"$0$",r"$\frac{F\ell}{EA} + \frac{n \ell^2}{2E}$",r"$\frac{F\ell}{2EA} +\frac{3 n \ell^2}{8E}$"])

# Draw background grid
ax[0].grid(True)

# Add a second subplot (you can customize this as needed)
# For demonstration, let's just add a simple plot in the second subplot
fem_sol_N=np.array([F+3/4*n*A*l,F+1/4*n*A*l])
ax[1].plot(x, F+n*A*(l-x), color='red', linewidth=2.0)
ax[1].plot([0,l/2],[fem_sol_N[0],fem_sol_N[0]],color='blue',linewidth=2.0)
ax[1].plot([l/2,l],[fem_sol_N[1],fem_sol_N[1]],color='blue',linewidth=2.0)
ax[1].set_xlabel(r"$x$")
ax[1].set_ylabel(r"$N(x)$")
ax[1].set_title("Schnittkraftverlauf")
ax[1].grid(True)
ax[1].set_xticks(tick_positions)
ax[1].set_xticklabels(tick_labels)

ax[1].set_yticks([1,F+n*l/2*A,F+n*l])
ax[1].set_yticklabels([r"$F$",r"$F+\frac{nA\ell}{2}$",r"$F+nA\ell$"])

# Create a common legend below both subplots
fig.legend(loc='lower center', bbox_to_anchor=(0.5, -0.0), ncol=2)

# Adjust layout to make room for the legend
#plt.tight_layout(rect=[0, 0.1, 1, 1])  # Adjust the rect to make space for the legend

# Show the plot
plt.show()

# Add legend.
# ax[0].legend([r"analytische Lösung", r"FEM"])
../../_images/176b6b3acfe78592c29b5be009b73e235cd8c000ef37a6aeb6a85d6e419913c3.png
from sympy import *
init_printing(use_unicode=True)
E,A,l,n = symbols('E A l n')
u_1, u_2, u_3 = symbols('u_1 u_2 u_3')
F, R = symbols('F R')
F_ext = Matrix([[F],[0],[R]])
F_v = (1/4)*n*l*A*Matrix([[1],[2],[1]])
F_g = F_ext+F_v
K=(2*E*A)/l*Matrix([[1,-1,0],[-1,2,-1],[0,-1,1]])
u=Matrix([[u_1],[u_2],[0]])
eq = F_ext+F_v - K*u
eq[0]
../../_images/560426c0d3df197dffa77059f36c13c287f3a9bacdb142a0fb4d5bbe4b041104.png
Kuu=K[0:2,0:2]
Kub=K[0:2,2]
Kbu=K[2,0:2]
Kbb=K[2,2]
F_g[0:2]
../../_images/0b16784cdfd5b684258c1a5e0723df372961e2379160ce876da9bdee56fe0d59.png
Kuu
../../_images/9154cf0a0710a4bf4854b3fd46f517859d0e1c4b6d7ad9dc58968c01107562e2.png
rhs=Matrix(F_g[0:2])-Kub*u[2]
rhs
../../_images/c4d9bc728640ecfd8bf2366691575049872f33c75544de19e4efaa6b69d2fc22.png
rhs=Matrix(F_g[0:2])-Kub*u[2]
display("rhs:",rhs)
uu = Kuu.LUsolve(rhs)
display("solution:",Kuu.LUsolve(rhs))
'rhs:'
../../_images/c4d9bc728640ecfd8bf2366691575049872f33c75544de19e4efaa6b69d2fc22.png
'solution:'
../../_images/82cbebb873944ef7dd02830693d99dfabeb32ac2247c1488fd7dd89759be5013.png
Kbu*uu-Matrix([F_v[2]])
../../_images/44a9584a1f3638edd5ecd19a92943fe03ac8a416eaa2a97d25afbfa8ac5a3518.png
Ke=E*A/l*Matrix([[1,-1],[-1,1]])
Fve=(1/4)*n*l*A*Matrix([[1],[1]])
ue=Matrix([[uu[1]],[uu[0]]])
Ke*ue-Fve
../../_images/5f5ec184a66a47c616c7960a7360b5deb66c8520d6fdeac1ee18192b2b9769e5.png
eps=2*E*A/l * Matrix([[1, -1]])
N1=-eps*Matrix([[0],[uu[1]]])
N2=-eps*Matrix([[uu[1]],[uu[0]]])
display("N1",N1)
display("N2",N2)
'N1'
../../_images/e65bfb1d8eea8ba70ab708147881ce5d9bbb752e0968b0b23253185b12106d52.png
'N2'
../../_images/c5d97ab60ac083e18494a1097940528574475efbd5ef570ce35c930ead4cddb6.png
uu
../../_images/82cbebb873944ef7dd02830693d99dfabeb32ac2247c1488fd7dd89759be5013.png

Stab plot#

import numpy as np
import matplotlib.pyplot as plt

# Set up the plot configurations.
plt.rcParams['lines.linewidth'] = 2.0
plt.rcParams['lines.color'] = 'black'
plt.rcParams['legend.frameon'] = True
plt.rcParams['figure.figsize'] = (8, 6)
plt.rcParams['font.family'] = 'serif'
plt.rcParams['legend.fontsize'] = 15
plt.rcParams['font.size'] = 15
plt.rcParams['axes.spines.right'] = False
plt.rcParams['axes.spines.top'] = False
plt.rcParams['axes.spines.left'] = True
plt.rcParams['axes.spines.bottom'] = True
plt.rcParams['axes.axisbelow'] = True
plt.rcParams['grid.color'] = 'grey'
plt.rcParams['grid.alpha'] = 0.6
plt.rcParams['grid.linestyle'] = '--'

sizex = 7.5
sizey = 5

# Create a figure and axis.
fig, ax = plt.subplots(1, 2,figsize=(2*sizex,sizey))

F=1
E=1
A=1
l=20
n=1

x = np.linspace(0, l, 100)


# Plotting the functions.
ax[0].plot(x, (F/(E*A)+n*l/E)*x-(n/(2*E))*x**2,color='red',linewidth=2.0,label="analytische Lösung")

fem_sol=np.array([
 [0,0],
 [l/2,F*l/(2*E*A)+3*n*l**2/(8*E)],
 [l,F*l/(E*A)+n*l**2 / (2*E)]   
])

#ax[0].scatter(fem_sol[:,0],fem_sol[:,1],color="blue",marker="o",s=100,label="FEM")
ax[0].plot(fem_sol[:,0],fem_sol[:,1],color="blue",marker="o",markersize=5,label="FEM")
# Set labels for the axis
ax[0].set_xlabel(r"$x$")
ax[0].set_ylabel(r"$u(x)$")
ax[0].set_title("Verschiebungsverlauf")

# Set custom ticks for the x-axis with LaTeX symbols
tick_positions = [0, l/2, l]
tick_labels = [r"$0$", r"$\frac{\ell}{2}$", r"$\ell$"]
ax[0].set_xticks(tick_positions)
ax[0].set_xticklabels(tick_labels)
ax[0].set_yticks([0,F*l/(2*E*A)+3*n*l**2/(8*E),F*l/(E*A)+n*l**2 / (2*E)])
ax[0].set_yticklabels([r"$0$",r"$\frac{F\ell}{EA} + \frac{n \ell^2}{2E}$",r"$\frac{F\ell}{2EA} +\frac{3 n \ell^2}{8E}$"])

# Draw background grid
ax[0].grid(True)

# Add a second subplot (you can customize this as needed)
# For demonstration, let's just add a simple plot in the second subplot
fem_sol_N=np.array([F+3/4*n*A*l,F+1/4*n*A*l])
ax[1].plot(x, F+n*A*(l-x), color='red', linewidth=2.0)
ax[1].plot([0,l/2],[fem_sol_N[0],fem_sol_N[0]],color='blue',linewidth=2.0)
ax[1].plot([l/2,l],[fem_sol_N[1],fem_sol_N[1]],color='blue',linewidth=2.0)
ax[1].set_xlabel(r"$x$")
ax[1].set_ylabel(r"$N(x)$")
ax[1].set_title("Schnittkraftverlauf")
ax[1].grid(True)
ax[1].set_xticks(tick_positions)
ax[1].set_xticklabels(tick_labels)

ax[1].set_yticks([1,F+n*l/2*A,F+n*l])
ax[1].set_yticklabels([r"$F$",r"$F+\frac{nA\ell}{2}$",r"$F+nA\ell$"])

# Create a common legend below both subplots
fig.legend(loc='lower center', bbox_to_anchor=(0.5, -0.08), ncol=2)

# Adjust layout to make room for the legend
#plt.tight_layout(rect=[0, 0.1, 1, 1])  # Adjust the rect to make space for the legend

# Show the plot
plt.show()

# Add legend.
# ax[0].legend([r"analytische Lösung", r"FEM"])
../../_images/c132fbcde5a0a9cd010188f575e66ae1f2b7a6eb1dbd7b12723ec7d3079eaa54.png

Balken Plot#

Shape Functions#

import numpy as np
import matplotlib.pyplot as plt
import scienceplots

# Set up the plot configurations.
plt.rcParams['lines.linewidth'] = 2.0;
plt.rcParams['lines.color'] = 'black';
plt.rcParams['legend.frameon'] = True;
plt.rcParams['figure.figsize'] = (8, 6);
plt.rcParams['font.family'] = 'serif';
plt.rcParams['legend.fontsize'] = 15;
plt.rcParams['font.size'] = 15;
plt.rcParams['axes.spines.right'] = False;
plt.rcParams['axes.spines.top'] = False;
plt.rcParams['axes.spines.left'] = True;
plt.rcParams['axes.spines.bottom'] = True;
plt.rcParams['axes.axisbelow'] = True;
plt.rcParams['grid.color'] = 'grey';
plt.rcParams['grid.alpha'] = 0.6;
plt.rcParams['grid.linestyle'] = '--';

plt.style.use(['science', 'grid']);


fig,ax = plt.subplots(1,1);

# Define the x values
x = np.linspace(0, 1, 100);
# Define the shape functions - Hermite polynomials 4th order
N1 = 1 - 3*x**2 + 2*x**3;
N2 = x - 2*x**2 + x**3;
N3 = 3*x**2 - 2*x**3;
N4 = -x**2 + x**3;
# Plot the shape functions
ax.plot(x, N1, label='$N_1$');
ax.plot(x, N2, label='$N_2$');
ax.plot(x, N3, label='$N_3$');
ax.plot(x, N4, label='$N_4$');
# Add labels and legend
ax.set_xlabel(r'$\eta$');
ax.set_ylabel(r'$N(\eta)$');
ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left');
# Show the plot
plt.show()
../../_images/701440f17342ebe2f3daa08e467b90d1d9404576f1b03a1dd899ecbb5af74619.png
import numpy as np
import plotly.graph_objects as go

# Define the x values
x = np.linspace(0, 1, 100)

# Define the shape functions - Hermite polynomials 4th order
N1 = 1 - 3*x**2 + 2*x**3
N2 = x - 2*x**2 + x**3
N3 = 3*x**2 - 2*x**3
N4 = -x**2 + x**3

# Create the figure
fig = go.Figure()

# Add traces for each shape function
fig.add_trace(go.Scatter(x=x, y=N1, mode='lines', name='N_1'))
fig.add_trace(go.Scatter(x=x, y=N2, mode='lines', name='N_2'))
fig.add_trace(go.Scatter(x=x, y=N3, mode='lines', name='N_3'))
fig.add_trace(go.Scatter(x=x, y=N4, mode='lines', name='N_4'))

# Update layout with LaTeX math
fig.update_layout(
    title=r'$\text{Hermite Polynomials 4th Order}$',
    xaxis_title=r'$\eta$',
    yaxis_title=r'$N(\eta)$',
    legend=dict(x=1.05, y=1),
    font=dict(family="Serif", size=15),
    template='plotly_white',
    xaxis=dict(showgrid=True, gridcolor='grey', gridwidth=0.6, griddash='dash'),
    yaxis=dict(showgrid=True, gridcolor='grey', gridwidth=0.6, griddash='dash')
)

# Show the plot
fig.show()
import sympy as sp
from sympy import Matrix, latex

# Define the variables
xi,ell = sp.symbols('xi, ell')

# Define the Hermite polynomials
N1 = 1 - 3*xi**2 + 2*xi**3
N2 = (xi - 2*xi**2 + xi**3)*ell
N3 = 3*xi**2 - 2*xi**3
N4 = (-xi**2 + xi**3)*ell

N = Matrix([[N1], [N2], [N3], [N4]])
dNdxi = sp.diff(N, xi)
dNdxidxi = sp.diff(dNdxi, xi)

M=dNdxidxi*dNdxidxi.T
M
../../_images/d57c4b7f7ed46d3f2cd73d1c0136e95905b517c6e7611060cac5148067b40ca8.png
print(latex(M))
\left[\begin{matrix}\left(12 \xi - 6\right)^{2} & \ell \left(6 \xi - 4\right) \left(12 \xi - 6\right) & \left(6 - 12 \xi\right) \left(12 \xi - 6\right) & \ell \left(6 \xi - 2\right) \left(12 \xi - 6\right)\\\ell \left(6 \xi - 4\right) \left(12 \xi - 6\right) & \ell^{2} \left(6 \xi - 4\right)^{2} & \ell \left(6 - 12 \xi\right) \left(6 \xi - 4\right) & \ell^{2} \left(6 \xi - 4\right) \left(6 \xi - 2\right)\\\left(6 - 12 \xi\right) \left(12 \xi - 6\right) & \ell \left(6 - 12 \xi\right) \left(6 \xi - 4\right) & \left(6 - 12 \xi\right)^{2} & \ell \left(6 - 12 \xi\right) \left(6 \xi - 2\right)\\\ell \left(6 \xi - 2\right) \left(12 \xi - 6\right) & \ell^{2} \left(6 \xi - 4\right) \left(6 \xi - 2\right) & \ell \left(6 - 12 \xi\right) \left(6 \xi - 2\right) & \ell^{2} \left(6 \xi - 2\right)^{2}\end{matrix}\right]
print(latex(N))
\left[\begin{matrix}2 \xi^{3} - 3 \xi^{2} + 1\\\ell \left(\xi^{3} - 2 \xi^{2} + \xi\right)\\- 2 \xi^{3} + 3 \xi^{2}\\\ell \left(\xi^{3} - \xi^{2}\right)\end{matrix}\right]
K=sp.integrate(M, (xi, 0, 1))
K
../../_images/332d0267287612f1ae58b718aa253b0bc5d3f0d883bdda048e3f6199cc4f78e6.png
print(latex(K))
\left[\begin{matrix}12 & 6 \ell & -12 & 6 \ell\\6 \ell & 4 \ell^{2} & - 6 \ell & 2 \ell^{2}\\-12 & - 6 \ell & 12 & - 6 \ell\\6 \ell & 2 \ell^{2} & - 6 \ell & 4 \ell^{2}\end{matrix}\right]
Nq = Matrix([[1-xi], [xi]])
q1,q2 = sp.symbols('q1, q2')
q = Matrix([[q1], [q2]])

Fq1 = N*(Nq.T*q)
Fq1
../../_images/94b24e924a5a0078788e4404387e2f55965d1d676526bf086b02d781a388f47d.png
print(latex(Fq1))
\left[\begin{matrix}\left(q_{1} \left(1 - \xi\right) + q_{2} \xi\right) \left(2 \xi^{3} - 3 \xi^{2} + 1\right)\\\ell \left(q_{1} \left(1 - \xi\right) + q_{2} \xi\right) \left(\xi^{3} - 2 \xi^{2} + \xi\right)\\\left(- 2 \xi^{3} + 3 \xi^{2}\right) \left(q_{1} \left(1 - \xi\right) + q_{2} \xi\right)\\\ell \left(\xi^{3} - \xi^{2}\right) \left(q_{1} \left(1 - \xi\right) + q_{2} \xi\right)\end{matrix}\right]
Fqi = sp.integrate(Fq1, (xi, 0, 1))
Fqi
../../_images/9b62cb0f72aa13ecd4f37fee6fd1aefea6c3683b7c7d6e170e40c88a31ddb1f2.png
print(latex(Fqi))
\left[\begin{matrix}\frac{7 q_{1}}{20} + \frac{3 q_{2}}{20}\\\frac{\ell q_{1}}{20} + \frac{\ell q_{2}}{30}\\\frac{3 q_{1}}{20} + \frac{7 q_{2}}{20}\\- \frac{\ell q_{1}}{30} - \frac{\ell q_{2}}{20}\end{matrix}\right]

Balken unter konstantanter Streckenlast#

from sympy import Rational
from sympy import lambdify
from sympy import Eq

x,EI,ell,q0,C_1,C_2,C_3,C_4 = sp.symbols('x,EI,ell,q0,C_1,C_2,C_3,C_4')

w = 1/EI*(Rational(1,24)*q0*x**4 +Rational(1,6)*C_1*x**3+Rational(1,2)*C_2*x**2+C_3*x+C_4)
wfun = lambdify(x, w)
w 
../../_images/44f4783ddf297a79cae0b174e59068d583e572d300e6f6624edfc0c1eae3735d.png
eq1 = Eq(w.subs(x,0),0)
eq2 = Eq(w.subs(x,ell),0)
eq3 = Eq(-w.diff(x,2).subs(x,0),0)
eq4 = Eq(-w.diff(x,2).subs(x,ell),0)
solution = sp.solve((eq1,eq2,eq3,eq4),(C_1,C_2,C_3,C_4))
solution
../../_images/5661e22268fa2386b7c28e87ce6b148b6e54a0105a1888082efe5146fb69988b.png
w.subs(solution)*EI
../../_images/705c5937e111f41b45805db17c9367c34ccc50ee190183ef6579279c4a127273.png
print(latex(w.subs(solution)*EI))
\frac{\ell^{3} q_{0} x}{24} - \frac{\ell q_{0} x^{3}}{12} + \frac{q_{0} x^{4}}{24}
print(latex(-w.subs(solution).diff(x,3)*EI))
q_{0} \left(\frac{\ell}{2} - x\right)
-w.subs(solution).diff(x,3)*EI
../../_images/3dc03f7c99758ed761290e0f12f82aabda95ab5d4be1a0fe577a4126686f1f2e.png
wp
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[33], line 1
----> 1 wp

NameError: name 'wp' is not defined
from sympy.plotting import plot, PlotGrid

wp = w.subs(solution).subs({"ell":1,"q0":1,"EI":1})
M = -wp.diff(x,2)
Q = -wp.diff(x,3)

p1= plot(-wp,(x,0,1),title='Durchbiegung $EI w$',show=False,xlabel=r'$\frac{x}{\ell}$',ylabel='$EI w$',line_color='black')
p2= plot(M,(x,0,1),title='Biegemoment $M$',show=False,xlabel=r'$\frac{x}{\ell}$',ylabel='$M$',line_color='red')
p3= plot(Q,(x,0,1),title='Querkraft $Q$',show=False,xlabel=r'$\frac{x}{\ell}$',ylabel='$Q$',line_color='blue')


PlotGrid(1,3,p1,p2,p3,size=(15,5))
../../_images/865a0f74f3c73f11465880442ab87fdde2b541def6b59c60a8272746b19a45c3.png
<sympy.plotting.plotgrid.PlotGrid at 0x7f04740a5de0>
wc = w.subs(solution)
eq1 = Eq(wc.diff(x,1),0)
sol = sp.solve(eq1,x)
sol
[ell/2, ell/2 + sqrt(3)*ell/2, -sqrt(3)*ell/2 + ell/2]
wc.subs({"x":ell/2})
\[\displaystyle \frac{5 \ell^{4} q_{0}}{384 EI}\]